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We investigate systems of identical bosons with the focus on two-body correlations and attractive 
finite-range potentials. We use a hyperspherical adiabatic method and apply a Faddeev type of 
decomposition of the wave function. We discuss the structure of a condensate as function of particle 
number and scattering length. We establish universal scaling relations for the critical effective 
radial potentials for distances where the average distance between particle pairs is larger than 
the interaction range. The correlations in the wave function restore the large distance mean-field 
behaviour with the correct two-body interaction. We discuss various processes limiting the stability 
of condensates. With correlations we confirm that macroscopic tunneling dominates when the trap 
length is about half of the particle number times the scattering length. 
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I. INTRODUCTION 

Condensation of a macroscopic number of bosons in 
the same quantum state was predicted many years ago 
[1]. Much later this was experimentally achieved in the 
laboratory for dilute systems of alkali gases [2-4]. The 
average properties of these gases are accounted for by 
the Gross-Pitaevskii equation [5, 6]. Exhaustive reviews 
of the theoretical developments after the experimental 
breakthrough can be found in [7, 8]. 

Degrees of freedom beyond the mean-field are crucial 
for the stability of the condensates, e.g., recombination 
into bound dimer and trimer cluster states [9-11]. The 
importance of such correlations is revealed in recent ex- 
periments [12, 13]. One example is the collapse of a Bose 
gas with large scattering length [13] where the lack of 
atoms in the condensate challenged the mean-field de- 
scription in terms of the time-dependent Gross-Pitaevskii 
equation [14, 15]. 

The mean- field description is valid for n|a s | 3 <C 1, 
where n is the density and a s is the two-body s-wave 
scattering length, when the particles on average are out- 
side the interaction volume of the order of scattering 
length, a s , to the third power [8]. The mean-field method 
neglects all correlations and thus breaks down at larger 
densities where correlations become important. Going 
beyond the mean-field is often very complicated as ex- 
emplified by Jastrow theory [16-18], which leads to high- 
dimensional equations. A later formulation is contained 
in the Faddeev- Yakubovskii equations [19, 20] where the 
wave function is expressed in terms of components descri- 
bing the asymptotic behaviour of all kinds of clusters. 

Comparison of different models is not always straight- 
forward, since different degrees of freedom are treated 
and the two-body interactions must be renormalized ac- 
cordingly. For example, using realistic potentials in self- 
consistent mean-field calculations leads to disastrous re- 
sults because the Hilbert space does not include correla- 
tions [21] as needed to describe both the short- and long- 
range asymptotic behaviour. In Gross-Pitaevskii calcu- 
lations the (5-function interaction is renormalized to give 



the correct scattering length in the Born approximation 
[8]. However, this substitution is only valid in the low- 
density limit. When correlations are included a different, 
and more realistic, interaction must also be used. Fur- 
thermore, only average properties can be described in the 
mean-field approximation. Thus comparisons of correla- 
tion dependent quantities are meaningless. 

Five years ago an interesting alternative study of a con- 
densate was formulated in terms of hyperspherical coor- 
dinates without any two-body correlations [22]. Using 
the same coordinate system a theoretical frame for de- 
scribing correlations was given soon after [23]. Detailed 
three-body calculations with zero total angular momen- 
tum were recently performed in the same framework [24] . 
Here the scattering length is varied and excited three- 
body states and any number of bound two-body molecu- 
lar states are allowed. The claim in [24] is that higher- 
lying Bose-Einstein condensed states in a trap of length 
b t do not collapse when N\a s \/b t > 0.5 as otherwise in- 
dicated by experiments [13]. The recombination takes 
place at distances several times the scattering length. 
They conjecture that the properties for N > 3 are quan- 
titatively similar to these three-body results. Another 
study in the same framework investigates the model- 
dependence of the three-body energy and finds that only 
the large scattering length enter [25]. They also indicate 
that the energy is insensitive to possible higher-order cor- 
relations for systems with many particles in a trap. 

In a further development using the adiabatic hyper- 
spherical expansion we formulated a method to describe 
two-body correlations in many-boson systems [26, 27]. 
This method is a novel attempt to describe correlated 
systems of low density. The formulation heavily relies on 
an additive set of components of the wave function as in 
the Faddeev decomposition but in contrast to the Jastrow 
multiplicative formulation. The numerical studies were 
limited to Bose-Einstein condensation for 20 particles. 

The purpose of the present paper is to extend the appli- 
cations to arbitrary scattering lengths and large particle 
numbers. We want to extract the general properties of 
the solutions especially for large scattering lengths where 
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mean-field computations are invalid. We obtain natu- 
rally self-bound many-body systems, even when the two- 
and three-body subsystems are unbound. The paper be- 
gins with a brief description of the hypersphcrical adia- 
batic expansion method in section II. Then in section III 
the details of the properties of the angular eigenvalue is 
discussed as the decisive ingredient in the radial poten- 
tial prociding the information about the two-body inter- 
action. In section IV we discuss the radial potential and 
the properties of the corresponding solutions. Finally in 
section V we discuss stability criteria for condensates ex- 
pressed in terms of various time scales and decay rates. 
Section VI contains the conclusions. 



II. HYPERSPHERICAL ADIABATIC METHOD 

We use the hyperspherical adiabatic expansion method 
with finite-range two-body interactions and simplifying 
assumptions about the wave function. We shall briefly 
describe the method and the assumptions. Details are 
given in [27]. 

The system of N identical particles of mass m is in the 
center of mass frame described by hyperspherical coordi- 
nates, i.e., one length, the hyperradius p, given by 



N N 



a) 



i<j 



and 3iV— 4 hyperangles [23, 28]. The ith single-particle 
coordinate is n, R is the center of mass coordinate, and 



nj = \ fj - n \ = \f2psh\a tj , 



(2) 



with ctij G [0, 7r/2] . The atoms are trapped in an external 
field approximated by a spherically symmetric harmonic 
oscillator potential of angular frequency u>: 



Kxt = £ Wrf = \rmj\p 2 + NR 2 ) • 



(3) 



i=l 



Without any two-body interaction between the partic- 
les the ground-state wave function is a Hartree product 
of Gaussian amplitudes: 



N 



total 



JJ e -</(2fc t 2 ) = e -p 2 /(2bi) e -NR'/(2b() ) (4) 



where the trap length is bt = y/h/ (mio). The second 
radial moments are (rf ) = 36 t 2 /2 and (R 2 ) = 3b 2 /(2N). 
For large N the average hyperradius therefore approaches 
the average mean-field radial coordinate times V~N ', see 
eq. (1). The hyperangles determine the relative orien- 
tations of the particles. 

With a two-body interaction term V(rij) the total 
Hamiltonian becomes 



N ~2 N 
i—1 i<j 



(5) 



It separates into a center of mass part (-ff c .m.), a radial 
part (Hp), and an angular part (/in) depending respec- 
tively on R, p, and f2 [28]: 
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(6) 

(7) 
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where T p and Tq are radial and angular kinetic energy 
operators. Then the center of mass motion always se- 
parates from the relative motion since the Vij -terms are 
independent of R. 

We remove the center of mass motion and study the 
Schrodinger equation for relative coordinates 



[H - H c . m .)V = . 



(10) 



The adiabatic hyperspherical expansion of the wave func- 
tion is 

oo 

*(p,n) = p-^ N -^' 2 Y,Up)^{p^), (ii) 

where $„ is an eigenfunction of the angular part of the 
Hamiltonian with an eigenvalue h 2 \ u (p) / [2mp 2 ): 



ftn*„(p,n) = A„(p)$„(p,ft) 



(12) 



Then eq. (10) leads to a set of coupled radial equations. 
Neglecting couplings between the different ^-channels 
yields the radial eigenvalue equation: 



2mU v {p) _ X v (37V - 4)(3A^ - 6) 



(13) 
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where E v is the energy and the adiabatic potential U v 
acts as an effective mean-field potential as a function of 
the hyperradius. This potential consists of three terms, 
i.e., the external field, the generalized centrifugal barrier, 
and the angular average of the interactions and kinetic 
energies. The neglected non-diagonal terms are typically 
about 1% of the diagonal terms for attractive Gaussian 
potentials. 

We have so far no restriction on the many-body wave 
function, but include in principle any structure of the 
system. To choose a convenient form we follow the philo- 
sophy in the Faddeev-Yakubovskil formulations [19, 20], 
i.e., the additive decomposition of the wave function re- 
flects explicitly the possible asymptotic large-distance be- 
haviour of cluster subsystems. We expect that two-body 
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correlations are most important and we select the cor- 
responding terms in the decomposition. Higher-order 
correlations are then essentially neglected. This proce- 
dure assumes a very different starting point compared to 
the Jastrow factorization into products of two-body wave 
functions [16, 29, 30]. The traditional Jastrow form is ex- 
pected to be more efficient for large densities while our 
method is well suited for the low densities encountered 
for Bose-Einstcin condensates. 

Emphasizing two-body correlations we therefore de- 
compose the angular wave function $ in the symmetric 
Faddeev components <p 

N N 

= XJ ^ ~ X <t>(P, r tj ) , (15) 

i<j i<j 

where the last approximation assumes that only rela- 
tive s-waves between each pair of particles contribute. 
Then the coordinate dependence reduces to the distance 
Tij = \[2p sin oiij. Neglecting higher-order partial waves 
is justified when the large- distance properties are deci- 
sive. The capability of this assumption for large scatte- 
ring length has been demonstrated for N = 3 by descri- 
bing the intricate Efimov effect [31, 32]. 

The angular eigenvalue equation (12) can by a varia- 
tional technique be rewritten as a second order integro- 
differential equation in the variable a\2 [28]. For atomic 
condensates the interaction range is very short compared 
to the spatial extension of the iV-body system. Then 
this equation simplifies even further to contain at most 
one-dimensional integrals. The validity of our approxi- 
mations only relies on the small range of the potential, 
whereas the scattering length can be as large as desired. 

We shall use the finite-range Gaussian potential 
y{ r ) — Vb exp(— r 2 /6 2 ). Thus we have either overall 
attractive or overall repulsive potentials depending on 
the sign of the strength Vq. It is convenient to mea- 
sure the strength of the interaction in units of the Born- 
approximation of the scattering length 

^ 4^ J d fkl V ™ = ' (16) 

where the last expression is for the Gaussian potential. 
We use the sign convention that the scattering length 
a s > for a purely repulsive potential, such that a s ~ 
ae for |ae|/6 <C 1. Thus a s > for purely repulsive 
potentials while purely attractive potential can lead to 
any, positive or negative, value of a s depending on Vq and 
b. In appendix A is collected the connections between the 
Gaussian strength measured in a^/b and the scattering 
length a s /b for the cases applied in this work. In most of 
the numerical work we have |ob|/6 close to unity. 

III. ANGULAR POTENTIALS 

The key quantity in the radial equation (13) is the 
angular eigenvalue A obtained from eq. (12). This eigen- 
value depends on the number of particles, on the size of 



the system through the hyperradius, and on the two-body 
potential through the scattering length. The behaviour 
of A is decisive for the effective potential in eq. (14) 
which in turn determines the properties of the solutions 
to eq. (13). We shall therefore first study the dependence 
of A on the parameters in the model. We use the method 
described in [27]. The two-body interaction is a sim- 
ple Gaussian either purely attractive or purely repulsive. 
This finite-range interaction never produces the collapse 
at short distance arising from an attractive <5-force [33]. 
Thus we can as well use attractive potentials with one or 
more bound states. 



A. General eigenvalue behaviour 

The angular eigenvalue spectrum coincides with the 
free spectrum (without interaction) at both small and 
large hyperradii; for p = because all interactions are 
multiplied by p 2 , see eqs. (6) and (9), and at p = oo be- 
cause the short-range interaction has no effect at large 
distances. Thus, perturbation theory for small p for a 
Gaussian potential shows that the eigenvalues change 
from their hyperspherical values A„(0) = 2^(2^+3^ — 5), 
v = 0, 1, . . ., as 

K{p)-K{Q) = ^N(N-l)p 2 . (17) 

If the two-body potential is attractive, but too weak 
to support any bound state, the eigenvalues reach a mi- 
nimum as function of p and then return to one of the 
finite hyperspherical values. For more attractive poten- 
tials there is a one-to-one correspondence between one 
given two-body bound state of energy < and one 
eigenvalue A diverging with p as A = 2mE^ 2 ' p 2 / h 2 . The 
corresponding structure describes, appropriately sym- 
metrized, one pair of particles in that bound state and all 
others far apart from the pair and from each other. In ad- 
dition to this finite number of such negative eigenvalues 
the hyperspherical spectrum emerges at large distances. 

To illustrate we show in fig. 1 a number of possible 
angular eigenvalues A as functions of hyperradius for dif- 
ferent potentials. The entirely positive (solid) curve cor- 
responds to a repulsive Gaussian. The diverging (dotted 
and thick dot-dashed) curves correspond to potentials 
with one bound two-body state. For our purpose the 
curves approaching zero for large p (dashed and thin dot- 
dashed curves) are the most interesting, since they are 
crucial for the later description of the condensate. This is 
true even when the potential has lower-lying bound states 
corresponding to diverging A (thick dot-dashed curve). 

The convergence of A as p — > is due to the finite 
range of the potential and the behaviour depends on the 
interaction range b. The deep minima in fig. 1 at small to 
intermediate distances depend strongly on both the num- 
ber of particles and the strength of the attraction. They 
are substantially deeper than reported in [26, 27] where 
one term inadvertently was used with the wrong sign in 
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FIG. 1: Angular eigenvalues A„ (numbered with inceasing 
v as v — 0,1) as functions of hyperradius divided by inter- 
action range, p/b, for N = 100, for different scattering lengths 
a s /b and numbers of bound two-body states A/"b indicated as 
\ v (a s /b, A/"b) on the figure. 



the numerical examples. This error is only significant 
for small and intermediate p. Increasing the strength of 
the attraction always leads to larger negative values of A. 
However, at some point one more bound two-body state 
reveals its presence by changing convergence to zero into 
a parabolic divergence with p. 

For distances much larger than the range of the po- 
tential the eigenvalues could as well be computed from 
a zero-range interaction, i.e., 4Trh 2 a s S(f)/m. The hyper- 
harmonic angular wave function should then be appro- 
priate for $ and the eigenvalue A obtained as the corre- 
sponding expectation value. The lowest hyperharmonic 
is a constant independent of angles and the result is [22] 



<2 TQ 
2 V 7T p 



N(N - 1) 



p 



(18) 



This zero-range result is inversely proportional to p for 
all hyperradii and consequently with a non-physical di- 
vergence when p — ► 0. The only length scale arises from 
the strength of the 5-function. In mean-field calculations 
this strength is chosen to reproduce the correct scatte- 
ring length as in the Born approximation [8, 21]. To 
reach this limit with a Gaussian potential then requires 
that the (5-function is approached while ob is maintained 
equal to the desired value of a s . 

This artificial construction is due to the lack of cor- 
relations in mean-field computations where the effective 
interaction is adjusted to the available Hilbert space. 



We use finite-range Gaussian potentials and include two- 
body correlations. Then we expect the large-distance 
asymptotic behaviour to be described by eq. (18) with 
the correct scattering length. This tests the efficiency of 
the simplified structure of the wave function in eq. (15). 
Mathematically this should result from the structure of 
the second order integro-diffcrcntial angular eigenvalue 
equation [27, 28]. 

Numerically we investigate the asymptotic behaviour 
of A in this context by comparing to the zero-range re- 
sult A<5 in fig. 2. The convergence to the limiting value 
is fastest for the smallest value of \a s \ (dashed and solid 
curve) reflecting that the correlations arising for large 
scattering lengths (dotted line) cannot be accounted for 
by the zero-range result. This is well understood for 
three particles where the Efimov effect (very large a s ) 
extends correlations in hyperradius to distances around 
four times the average scattering length [31, 32]. These 
effects are not present in the zero-range expectation value 
contained in A,5. When p exceeds \a s \ by a sufficiently 
large amount the Efimov effect disappears in A and Xs is 
approached. 
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FIG. 2: Same as figure 1, but the angular potential is shown 
in units of the zero-range result in eq. (18) as obtained in [22]. 

A stronger attraction corresponding to one two-body 
bound state produces one diverging eigenvalue (figure 1) 
while the second eigenvalue converges towards Xs (dot- 
dashed curve). In fact Ai(— 10, 1) almost coincides with 
the lowest eigenvalue Ao(— 10,0) for the same scattering 
length but for a potential without bound two-body states 
(dotted curve). 

The numerical deviations from Xs at large distance is 
in all cases less than 10%. The asymptotic behaviour is 
very smooth but still originating in systematic numerical 
inaccuracies. 

These results demonstrate that the scattering length 
entirely determines the asymptotic behaviour of the 
potentials. The radial shape of the two-body poten- 
tial could be Gaussian, square-well, Woods-Saxon, or 
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Yukawa, still the same a s would produce the same an- 
gular eigenvalue at sufficiently large distance. 



B. iV-dependence 

The angular eigenvalues increase rapidly with TV as 
seen already from the 7V 7 / 2 -dependence in \$. The ma- 
jor variation in magnitude is then accounted for by using 
this large-distance zero-range result as the scaling unit. 
We show in fig. 3 a series of calculations for the same 
two-body interaction for different numbers of atoms. All 
curves are similar with a systematic increase in the char- 
acteristic hyperradius p a where they bend over and ap- 
proach the zero-range result. We then numerically deter- 
mine this characteristic length p a to be proportional to 
the scattering length and a particular power 7/6 of N, 
i.e., 

p a (N) = \a s \N 7 / 6 . (19) 




10 3 10 4 10 5 10 6 10 7 10 1 

p/b 



FIG. 3: The lowest angular eigenvalue in units of \g as a 
function of hyperradius for a s /b — —401 for four different 
numbers of particles N = 10 2 , 10 3 , 10 4 , 10 5 . 

The quality of this scaling is illustrated in fig. 4, where 
all curves essentially coincide for distances smaller than 
p a - At larger hyperradii the zero-range result of +1 
should be obtained. However, systematic deviations from 
a common curve is apparent. For each N one smooth 
curve is followed at small and intermediate distances im- 
plying that the numerical inaccuracies here are system- 
atic until random fluctuations set in at large p. 

The smooth numerical curves can be rather well repro- 
duced by the function 

\(-\N,p) = \ 5 (N,p)-g(-\p/p a ), (20) 
g<-\x) = 9oo (l-e-*/*-)(l + ^) , (21) 
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FIG. 4: The same as fig. 3, but with p in units of p a . The 
points following the intermediate curve (goo — 0.8) are ob- 
tained with many integration points and the points along the 
lower curve (goo = 0.4) are obtained with fewer points. The 
curve for goo — 1.0 is the expected correct asymptotic be- 
haviour from eq. (20). The points for N = 3 are calculated 
with the zero-range model from [34]. 

where goo — 1 in accurate calculations. The exponential 
term is introduced to reproduce the rather fast approach 
to the asymptotic value as seen in fig. 4. The behaviour 
at smaller distance, depending on the range of the inter- 
action, is simulated by the Xb-term. The extreme limit 
of p — > is attempted reproduced on the function in 
eq. (21). 

TABLE I: Numerical values of goo, x a , and Xb for four scat- 
tering lengths. 



a s /b 


-5.98 


-401 


-799 


-4212 


goo 


0.99 


0.80 


0.65 


0.30 




1.06 


0.74 


0.59 


0.28 


Qoo f %a 


0.93 


1.081 


1.099 


1.077 


Xb 


0.15 


2.3 • 10~ 3 


1.15 • 10~ 3 


2.2 • 10~ 4 


Xb/(b/\a s \) 


0.92 


0.922 


0.919 


0.927 



The two groups of computations in fig. 4 are reason- 
ably well reproduced by the parameter sets x a ~ 0.74, 
x ~ 2.3 • 10~ 3 , and g^o ~ 0.8 or g^ ~ 0.4. These pa- 
rameters may depend on the scattering length, and we 
therefore repeated the computation for various a s . The 
best choice of parameters are shown in table I. We no- 
tice that goo and x a both are of order unity, and that the 
fraction gco/x a is almost constant, except for the small- 
est scattering length. The parameter x b , introduced to 
account for the finite interaction range, is almost equal 
to 6/|a s |. At large hyperradii, where p 3> p a , A^ _ ) ap- 
proaches gooXs- The rather accurate results for N = 100 
displayed in fig. 2 confirm that g^ ~ 1 by deviating less 
than 10% from \$ at large hyperradii. 

The angular eigenvalue is given by g(~\x) ~ goox/x a 
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for xi, <C pj p a -C x a . Numerical calculations in this inter- 
mediate region of hypcrradii therefore rather accurately 
determines the fraction goc/x a ~ 1.08 as given in table I. 
With . 9oc = 1 this implies that x a ^ 1/1.08 ~ 0.92. The 
parameters of g(~\x) in cq. (21) are then given by 



1, 



0.92 



x b 



0.92— . (22) 

\&8 



We can compare with the rigorous result for N — 3 [32] 
where the angular eigenvalue at large-distance coincide 
with Xs in eq. (18), i.e., 



\ 5 (N = 3,p) = 



48a s 



(23) 



V2wp ' 

Thus also for N = 3 the universal function </~) asymp- 
totically approaches 1 for all scattering lengths. More 
accurate results for N = 3 have been calculated with 
the zero-range model from [34] and are shown in rescaled 
form in fig. 4. The behaviour is similar to the behaviour 
of the eigenvalue for the iV-body systems, which confirms 
the schematic model. 

The accuracy of the parametrization in eq. (21) is seen 
in figs. 5a-d, where the angular eigenvalues are shown in 
units of A^~) with the individual set of parameters from 
table I. Good agreement is found for p/p a > x 0l ex- 
cept at large hyperradii where the numerical inaccuracy 
increases with increasing scattering length. Fortunately, 
the large-distance behaviour is known from analytic con- 
siderations and we do not need to rely on numerical com- 
putations at these distances. The remaining deviations 
occur at small hyperradii for p/p a < x h or equivalently 
for p < N 7 ^ 6 b, where the result depends on the radial 
shape of the two-body interaction. 

C. iV-dependence with bound two-body states 

In the presence of a bound two-body state of energy 

£(2) 

one angular eigenvalue eventually diverges at large 
hyperradii as [35] 



A(2)(/0) = W^ (2) 



£ (2) < . 



(24) 



In the limit of weak binding, or for numerically large 
scattering lengths, the energy of the two-body bound or 
virtual state is given by 



h 2 



mat 



(25) 



where c approaches unity for large scattering lengths. 

We now parametrize the angular eigenvalue by an ex- 
pression similar to eqs. (20) and (21). The effect of 
the bound two-body state is only expected to show up 
at large distances where the behaviour corresponds to 
eq. (25). The small and intermediate distances resemble 
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FIG. 5: The lowest angular eigenvalue Ao in units of A'~', 
eqs. (20) and (21) and table I, as functions of the hyperradius 
in units of p a , eq. (19). The scattering lengths are given by 
a) a s /b = -401, b) a s /b = -799, c) a s /b = -4212, and d) 
a s /b = —5.98. The different iV-values are as indicated. 



the behaviour when no bound state is present. Therefore 
we arrive at the parametrization 



X^(N,p) = \ s (N,p) g( + \p/ Pa ), 



(26) 
(27) 



with the notation and estimates from eq. (22). 

We compare in fig. 6 the parametrization in eqs. (26) 
and (27) with the computed angular eigenvalues for a 
potential with one bound two-body state. For the large 
scattering length (a s /b = 100) in fig. 6a one smooth curve 
applies for all the particle numbers; numerical inaccura- 
cies set in at larger hyperradii, which is most obvious for 
the largest particle numbers. This smooth curve is in 
a large interval of hyperradii at most deviating by 20% 
from the parametrized form, and even less than 10% at 
large hyperradii, before the numerical instability sets in. 

For smaller a s (a s /b — +10) the deviation at large 
hyperradii is less than 1%. The deviation at intermediate 
distances would decrease by inclusion of a linear term in 
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eq. (27). The smooth curve at small hyperradii is outside 
the range of validity of the parametrization, i.e., within 
the range of the two-body potential and then depending 
on details of the interaction. 
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FIG. 6: a) The lowest angular eigenvalue Ao in units of X^ + \ 
eqs. (26) and (27), for a s /b = +100 and c = 1.02, when 
the potential holds one bound two-body state. The number 
of particles is indicated on the figure. The parameters are 
goo/x a = -1.09 and x b = 9.2 • 10~ 3 . b) The first excited 
angular eigenvalue Ai in units of As for a s /b = +10. 



inside and independent of the confining external field. 

This feature is absent in the mean-field description of 
Bose-Einstein condensation. For overall repulsive poten- 
tials corresponding to positive scattering lengths no at- 
tractive part is possible. For attractive potentials either 
a zero-range potential would produce a collapsed wave 
function and a finite-range potential would not give re- 
pulsion at large distance. 



D. The properties of A^' 

The functions X^' coincide when p <C p a and depends 
only on the absolute magnitude of the scattering length 
\a s \. For /)>/)„ the functions differ qualitatively, i.e., 
A(~) approaches zero as X$ while A^+) diverges as — p 2 . 

At intermediate hyperradii, xt, <p/p a «i a , when 



6< 



9 



N 7/6 



(28) 



the angular eigenvalue approaches a constant value 



Aoc — Xs~ 



Pa%a 



^N 7 '\[^~ -1.597V 7 / 3 . (29) 



This numerical result is in agreement with the following 
derivation. 

The angular eigenvalue for large scattering length a s is 
independent of hyperradius p when p is large compared 
to the range b of the potential but small compared to 
a s . The plateau value Aqo can be estimated as the in- 
tersection between two curves at the point p a . The first 
curve is the parabollically decreasing A(p) corresponding 
to a bound two-body state, i.e., X(p) = 2mp 2 E^/h 2 = 
-2p 2 /a 2 s , where is given by eq. (25) with c = 1. The 
second curve is the increasing A^ (p) for an attractive po- 
tential (a s < 0), see eq. (18). Thus Xs(p a ) — X(p a ) gives 



The lowest eigenvalue Ao diverges at large hyperradius 
as described by eq. (26). If the two-body potential only 
has one bound state the second eigenvalue Ai is expected 
to approach zero at large distances as A,5. This pat- 
tern should be repeated with more than one bound two- 
body state, i.e., the first non-divergent angular eigenvalue 
would behave as A^ for large p. 

We therefore in fig. 6b compare the computed first ex- 
cited angular eigenvalue with Aa for different N. As in 
fig. 4 we obtain smooth and almost universal curves at 
large p, where the approach to unity sets in exponen- 
tially fast depending on N, but now much later when 



Pa - 



;^V 7/6 |a s | , 



(30) 



Ao^AO ~ --^/|A 7 / 3 ~-1.65A 7 / 3 , (31) 

which is very close to the numerical results in eqs. (19) 
and (29). 

The symbol Aoo is chosen for this constant, since the 
/^-region where A = A^ increases proportional to \a s \, see 
eq. (28), and thus extends to infinity for \a s \ = oo. With 
no bound two-body states (a s < 0) the lowest angular 
eigenvalue approaches zero at larger hyperradii, whereas 



p ~ 10 p a - Clearly a parametrization would also here be it diverges towards — oo as p 2 when a bound two-body 



possible. 

The large-distance asymptotic behaviour of Ai now 
corresponds to an effectively repulsive potential. How- 
ever, at small and intermediate hyperradii the potential is 
still effectively attractive (Ai < 0). This attractive region 
may support a self-bound system located at distances far 



state is present (a s > 0). On the threshold for a two- 
body bound state a s = ±oo and the angular eigenvalue 
therefore remains constant. 

In [26] Aoo — —5 A 2 was estimated by courageous ex- 
trapolation of calculations for A = 10, 20, 30 and the 
analytic result for A = 3. The much better estimate in 
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eq. (29) of the large- TV asymptotics of Aoo (N) increases 
with a slightly higher power of N but with a smaller pro- 
portionality factor. 

IV. RADIAL POTENTIALS AND SOLUTIONS 

The radial equation is the next step in the process of 
obtaining knowledge about the physical properties of the 
many-boson system. The angular potentials found in the 
previous section now enter the effective radial potential 
in eq. (14) and infer information about the interactions to 
quantities like energy, size, and structure of the system. 

A. Properties of the radial potential 

The radial potential in eq. (14) consists of three terms 
where the repulsive centrifugal barrier and the confin- 
ing external field both are positive. The interaction 
term \ v can be either repulsive or attractive depending 
on hyperradius and which eigenvalue we consider. The 
combination has structure depending on the interaction. 
For a purely vanishing or repulsive two-body potential 
we arrive at a simple behaviour qualitatively similar to 
the non-interacting (dashed) curve shown in fig. 7 for 
N = 100. All solutions are confined to the region be- 
tween the infinitely large potential walls at small and 
large hyperradii. 

For a moderately attractive two-body potential a dif- 
ferent structure already appears for the lowest angular 
potential (solid curve in fig. 7). The large- distance be- 
haviour is determined by the trap and is roughly as with- 
out interaction, but the barrier at intermediate distance 
is now finite both in height and width. The barrier height 
is small compared to the potential at both small and large 
hyperradii. At smaller hyperradii a rather deep and re- 
latively narrow minimum is present outside a hard core 
repulsion. The minimum occurs for N = 100 at about 
150 times the range of the interaction which corresponds 
to a mean distance (2(p 2 ) /iV) 1 / 2 between each pair of 
particles of about 15 times the interaction range b. 

With this potential we solve the diagonal radial equa- 
tion. The solutions can be divided into groups related 
to either the first or the second minimum. The lowest- 
lying of the first group of solutions have negative ener- 
gies. In the model they are truly bound states as they 
cannot decay into continuum states at large hyperradii 
[26]. Their properties are independent of the external 
trap which only has an influence at much larger distances. 
These self-bound TV-body states can decay into lower- 
lying states consisting of various bound cluster states, 
e.g., diatomic or triatomic clusters. The possibility of 
self-bound many-body systems, even though the two- and 
three-body sub-systems are unbound, is also discussed 
by Bulgac [36], who, however, considers the three-body 
interaction strength as a determining parameter for the 
properties of the self-bound many-boson system. 




FIG. 7: a) Radial potential Uo from eq. (14) corresponding to 
the lowest angular potentials for N = 100 and a s /b — —1.0. 
We model the experimentally studied systems [12] of Rb- 
atoms with oscillator frequency v — lo/(2tt) = 205 Hz and 
interaction range b = 10 a.u., thus yielding &t = \Jh/(mui) = 
14426. Also shown as horisontal lines are the negative energies 
Eo,n, n = 0, . . . , 58 in the lowest potential in the uncoupled 
radial equation, eq. (13). b) Detail at larger hyperradii. The 
energy of the first oscillator-like state (see text) is shown as a 
horisontal line close to zero. 

The group of states in the higher-lying minimum at 
larger distance all have positive energies. They are only 
stable due to the confining effect of the external trap 
potential. The lowest of these is interpreted as the state 
of the condensate and indicated by a horizontal line in 
fig. 7b. This second minimum almost coincides with the 
minimum of the radial potential arising without any two- 
body interaction. Thus the structure of the condensate is 
similar for both positive and negative scattering lengths 
arising from either attractive or repulsive interactions. 
However, an attraction produces in addition a series of 
lower-lying states at smaller hyperradii. 

Increasing N leaves semi-quantitatively the same fea- 
tures for pure repulsion, whereas an unchanged attrac- 
tion leads to decreasing barriers at intermediate hyper- 
radius and at some point this barrier vanishes altogether. 
At the same time the attractive minimum at smaller 
hyperradius becomes deeper. This in turn leads to an 
increasing number of bound states in this minimum as 
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function of N. 

As the scattering length increases, the barrier dis- 
appears and the effective potential inside the trap has 
the p _2 -behaviour characteristic for Efimov states, see 
eq. (14) with A = Aoo of eq. (29). The lowest-lying states 
are influenced by the details of the two-body interaction 
and without Efimov features. However, the higher-lying 
states, located for p- values obeying eq. (28), exhibit the 
Efimov scaling. They easily become very large and lo- 
cated far outside the minimum responsible for the bind- 
ing. Only a finite number of bound states is possible due 
to the confining external field. 

These states are many-body Efimov states arising 
when the two-body scattering length is large. This is also 
precisely the condition for the three-body Efimov states 
[31, 37]. Therefore the many-body Efimov states are em- 
bedded in the continua of dimer, trimer and higher-order 
cluster states. They could be artifacts of the model where 
only special degrees of freedom are treated. However, 
these states may also be distinguishable resonance struc- 
tures which are relatively stable because the particles are 
very far from each other and the couplings to the contin- 
uum states therefore are very weak. So far this remains 
an open question. 



B. Interaction energy 

The total energy of a state in the first minimum arc 
independent of the external field as these states are lo- 
cated at small distances. These states have no analogue 
in mean-field calculations. In contrast, total energies of 
the states in the second minimum are dominated by the 
contribution from the confining field and therefore are 
rather insensitive to anything else than this field and the 
corresponding harmonic oscillator quantum numbers. It 
is then much more informative to compare the interaction 
energies where the large external field contribution is re- 
moved. 

In fig. 8 is shown the interaction energy per particle as 
a function of the particle number for a relatively weak 
attraction corresponding to a small scattering length. 
The Gross-Pitacvskii solution exists and the related in- 
teraction energy is negative due to the attraction between 
the particles. A nearly linear behaviour is observed at 
small particle numbers, since each particle interacts with 
the N — 1 other particles. As N increases the mean- 
field attraction increases and the Gross-Pitaevskii solu- 
tion becomes unstable when N\a s \/bt > 0.58. This cor- 
responds to N — 1000 with the present parameters of 
\a s \/b t = 58-10- 5 . 

These mean-field interaction energies are in fig. 8 com- 
pared to the results obtained with the correlated wave 
functions from the present formulation. Only a few of the 
large number of bound states for each A^-value resemble 
the Gross-Pitaevskii solutions with a radius correspond- 
ing to the second minimum. For N = 20 the lowest six 
states are located in the first minimum. Their interaction 
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FIG. 8: Gross-Pitaevskii interaction energy as a function of 
N for bt/b — 1442 for two a s -values. The points are results 
from the present work for N = 20, 100, 900 and a s /b = -0.84 
(ae/6 = —0.5). Only states in the second minimum are dis- 
played. 



energies are large and negative outside the scale of fig. 8. 
The seventh state is located in the second minimum with 
an interaction energy very close to the Gross-Pitaevskii 
result, while the eighth has a positive interaction energy. 
This feature is repeated for increasing N, i.e., the low- 
est state located in the second minimum is similar to 
the mean-field result and the higher-lying states in this 
second minimum are less bound. When the mean-field 
solutions collapse, the correlated solutions remain stable 
due to the use of a finite-range potential. 

The correlated and mean-field interaction energies are 
remarkably similar when both exist. It may at first 
appear odd that the mean-field interaction energy is 
marginally lower than by use of the better suited form 
of the correlated wave function. The reason is that we 
compare the mean-field result for an effective interaction 
which has the correct scattering length in the Born ap- 
proximation while the correlated solution is obtained for 
an interaction with the correct scattering length. The 
mean-field interaction is more attractive to compensate 
for the limited mean-field Hilbert space. The more re- 
vealing comparison is to use the same interaction in both 
calculations. 

We can then compare results for the same a^/b = 
—0.5, i.e., a Gross-Pitaevskii calculation with a s /b = 
—0.5 and a Gaussian of a^/b = —0.5 corresponding to 
a s /b = —0.84. As seen in fig. 8 (dashed curve) now 
the mean-field energies are much smaller. However, it 
is remarkable that the correlated solution essentially re- 
produces the energy of the mean-field calculation where 
the interaction is renormalized to reproduce the correct 
energy, but with the wrong wave function. The impli- 
cation is that the correlated wave function is sufficient 
to describe the correct structure with the correct inter- 
action. The large-distance average properties are at best 
obtained in mean-field computations, but all features of 
correlations are absent by definition. 
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C. Definition and size of a condensate 

The lowest-lying positive energy solutions located in 
the second minimum have properties similar to the con- 
densates obtained in Gross-Pitaevskii calculations. The 
present formulation also provides lower-lying negative- 
energy states. It is therefore necessary to discuss how to 
distinguish a condensate state from other (perhaps very 
unstable) TV-body states. 

In mean-field treatments, with repulsive two-body po- 
tentials and confining trap potentials, the condensate is 
uniquely defined as a statistical mixture of single-particle 
states with the ground state dominating [8, 38]. A con- 
densate has on average many particles in the lowest 
single-particle state. This many-body state is only stable 
against total fragmentation due to the trap and as such 
first of all determined by the properties of the trap. Even 
with the trap the many-body state is still at best only 
approximately stationary due to the neglected degrees of 
freedom which allow energetically favored di-, tri-, and 
multi-atomic cluster states. This instability is also an ex- 
perimental fact seen by permanent loss of trapped atoms, 
e.g., in recombination processes [13]. 

Without any two-body interaction the properties of the 
many-body system is determined by the thick, dashed po- 
tential curve in fig. 7. Then we can easily identify the con- 
densate as a state in this potential where the dominating 
component for finite temperature is the ground state. In- 
cluding attractive two-body interactions (full curve) the 
deep minimum at small hyperradius is produced. Then 
the corresponding ground state, located in this minimum, 
has nothing to do with a condensate. The density is so 
high that couplings to other degrees of freedom would de- 
velop higher-order correlations and processes like three- 
body recombinations would quickly destroy the single- 
atom nature of the gas. This iV-body ground state does 
not show the signature of a Bose-Einstein condensate, 
where many particles occupy one single-particle level. 

The formulation in the present work does not use the 
concept of single-particle levels. Therefore we cannot talk 
about a statistical distribution of particles with the ma- 
jority in the lowest state. However, we can talk about 
a many-particle system described as a superposition of 
many-body eigenstates, where the lowest states are fa- 
vored in thermal equilibrium. To clarify we can think of 
a quantum state $ as a superposition of different eigen- 
states ^ n (p,n) in cq. (11) given by 

oo 

*(p,n) = ^ Cn *„(p,o) 

n=0 

oo oo 

= p-^ N -^' 2 E c « E um*»<j>i , (32) 

n=0 l/=0 

with the normalization J2 n l c «| 2 = 1- A condensate must 
be sufficiently large to exceed a certain minimum inter- 
particle distance, d c , below which the atoms are too close 
and recombine very fast. This distance depends on the 



scattering length and on the number of particles. There- 
fore, in our formulation the stationary states cannot be 
characterized as a condensate if (r?) <C d% for this wave 
function. 

We define one of the stationary states in this model 
as the "ideal condensate" state, i.e., the state of lowest 
energy with 

(r%) nc > d c , (33) 

characterized by n = n c . This state is dominated by the 
component in the lowest adabatic potential although not 
necessarily the states of lowest energy, which might have 
an average particle distance less than d c . The appro- 
priate of these excited states depends on the number of 
particles and on the scattering length. The ideal conden- 
sate is then characterized by one dominating component, 
with |c„ c | ~ 1 and |c n ^„ c I < 1. 

If d c is significantly smaller than the trap length 6 t > 
then the state of lowest energy located in the second mi- 
nimum can be identified as the condensate. This state 
is characterized by a radial wave function f(p) with the 
root mean square radius (p 2 ) approximately equal to the 
hyperradius at the second minimum of the lowest adia- 
batic potential Uo(p). 

To be specific we show in fig. 9 the root mean square 
interparticle distance given by (r?-)„ = 2(p 2 ) n /(N— 1) for 
the lowest excited states (labeled by n) in the potential 
of fig. 7. All states with n < 58 have negative energy 
and 206 < {{r 2 l;j ) n ) 1/2 < 1006, which implies that the 
particles are separated more than their interaction range. 
Whether these average distances allow qualification as 
condensates depends on the decay rate of these states. 
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FIG. 9: The root mean square distance for v — as a, function 
of the quantum number for N = 100, a 3 /b — — 1, bt/b = 1442. 

When n > 59 the energies are positive and the av- 
erage particle distance now suddenly exceeds 20006. In 
fact we now find (rfj) ~ 36^ which approximately is ob- 
tained in the limit of a non-interacting gas. These states 
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probably qualify as condensates. Their interaction ener- 
gies are in fig. 8 compared to the Gross-Pitaevskii values. 
The discontinuity at n — 58, 59 is due to the intermediate 
barrier. Decreasing and eventual removal of the barrier 
would smear out this abrupt change of size. Some of the 
negative-energy states could then extend very far out and 
in fact have sizes comparable to the trap length. This 
investigation could be repeated for the higher adiabatic 
potentials, still neglecting the couplings. The same pat- 
tern is obtained with fewer states of small interparticle 
distance. 



V. DECAY RATES 

The condensate is unstable due to the neglected cou- 
plings into other degrees of freedom. The condensate 
therefore has to be located at relatively large distances. 
The decisive radial potentials are sensitively depend- 
ing on the scattering length. In fig. 10 we illustrate 
the different behaviour by using the angular eigenvalues 
parametrized through eqs. (14), (20), (21), (26), and (27). 
In fig. 10a the scattering length is relatively small and a 
large barrier separates the outer minimum from the inner 
region. By increasing the scattering length the barrier de- 
creases first into a relatively flat region as in fig. 10b and 
then disappears completely as in fig. 10c when the trap 
length is exceeded. With these potentials we can now 
discuss various decay processes, i.e., three-body recom- 
bination into dimers, macroscopic tunneling through the 
barrier and macroscopic collapse after sudden removal of 
the barrier. 




p/b 



A. Three-body recombination 

Bound state dimers can be formed by a three-body 
process where the third particle ensures conservation of 
energy and momentum. The number of these three-body 
recombination (rec) events per unit volume and time can 
be estimated by the upper limit given in [9, 11]: 



67.9 



H\a s 



m 



(34) 



where n is the density of the gas. This expression can 
be converted into an estimate of the recombination rate 
for a given hyperradius p. With the volume V = N/n, 
the relation between density and mean distance l/n = 
47r(r2.) 3 / 2 /3, and (r?-) = 2(p 2 )/(N - 1) obtained from 
eq. (1), the total recombination rate becomes 



rcc 

~~h 



v rec V ~ 0.5 



h\a s \ 4 N 4 
mp 6 



(35) 



where the mean square average p is defined as p 2 = (p 2 ). 
In the spirit of the adiabatic hyperspherical expansion 
method we use p as a classical parameter. The recombi- 
nation rate increases rapidly with decreasing p, as indi- 
cated by the vertical arrows in fig. 10. 



FIG. 10: The radial potential from the schematic model for 
N = 100, b t /b = 10 4 and a) a s /b = -6, b) a s /b = -50, c) 
a s /b — > — co. The wave function is the lowest radial solution 
in the non-interacting case. The horizontal lines in parts a) 
and b) indicate an energy level (not to scale). 



The recombination time T roc is defined by N(t) = 
N(0) exp(— i/Trcc), where N(t) is the number of remain- 
ing atoms. We then obtain r roc /ft = —dN/dt = N/T ICC 
and 



Nh 

r 

J- roc 



2mp 6 mrfj 
% S | 4 7V 3 = 4H\a s \ 



(36) 



where we used fy = ^J2/Np. The final expression for 
T roc is independent of N. Since the condensate has to 
form in the external trap it is reasonable to define stabil- 
ity against recombination by T rcc 3> T c = 2it/uj, where 
T c , is the oscillator time. With 1/u) = mb 2 /h and eq. (36) 

we get stability when fy » -v/STrlasI 2 / 3 ^ 3 = d c , which 
provides d c introduced in section IV C. In units of b t we 
obtain 



— = </8tt(— 
bt \ b t 



2/3 



(37) 
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Thus for |a s |/&t < 1 also d c /bt < 1. The average par- |a s |/6t <C 1. Therefore close to this threshold we have 



ticle distance for a state located in the second mini- 



mum is of the order b t and therefore 



these states is larger than the critical stability length 
d c . These states then qualify as condensates. For 87 Rb 
atoms (a s ~ 100 a.u.) trapped in a field of v ~ 100 Hz 
we obtain T rec ~ 7 days. 



B. Macroscopic tunneling 

The second decay process is related to macro- 
scopic tunneling through the barrier, as indicated in 
fig. 10b. The present model provides stationary eigen- 
states (within the allowed Hilbert space) which by de- 
finition are time independent. Thus, strictly the states 
do not tunnel through the barrier, but an exponentially 
small tail extends to small hyperradii. All particles in this 
tail would immediately recombine into molecular clus- 
ters, because the density is very large in the inner region 
(both p and r,j are small) . The rate of this two-step de- 
cay, i.e., tunneling through the barrier and subsequent 
recombination, can be computed as the knocking rate 
multiplied by the transmission coefficient, which is a mea- 
sure of the ratio of the probabilities at the turning points 
inside and outside the barrier. The rate of recombina- 
tion due to macroscopic tunneling can then be estimated 
semi-classically as in [22] by 



Nv 
1 + e 2ff 
1 



/ 1 d 2 U{p) 




m dp 2 


Pmin 



J Pin 



(38) 
(39) 
(40) 



where the factor N is needed to give the total number 
of recombined particles. Here p m i n is the position of the 
second minimum of U and p ln and p ut are the classical 
turning points of the barrier. 

The barrier depends strongly on the combination 
iV|a a |/&t [22, 26]. When N\a s \/b t < 1 the barrier is 
large and the very small rate can be estimated through 
eqs. (38), (39), and (40). The WKB action integral is 



a ~ -TV In 
2 



N\a s 



(41) 



The barrier is absent when N\a s \/bt > 0.53. Close to, 
but before reaching, this threshold of stability the WKB- 
exponent can be approximated by 



a ~ 1.7N6 S 



1 



N\a s \/b t 
0.53 



(42) 



which is valid when S s is close to zero. 

The barrier is observed to vanish when N\a s \/b t ~ 
0.53 [12, 13], which due to the factor of TV implies that 



for a condensate in the second minimum that 



b t > 



for d c , i.e., the three-body recombination does not limit the 



stability. In the limit a <C 1 we get explicitly 



1 



7.07V 4 



< 1 



(43) 



implying that the macroscopic tunneling process domi- 
nates. With it « 1 we obtain that T tun /h = N/T tun — 
O.SNis, which for v ~ 100 Hz corresponds to a macro- 
scopic tunneling time of 10 ms. This is much faster than 
the three-body recombination time when the barrier is 
small (a < 1), i.e., T rec > T tun , see eq. (43). 

The three-body recombination rate is in fig. 11 shown 
as a function of hyperradius (solid curve) and com- 
pared with the macroscopic tunneling rate (dashed curve) 
where all particles in the condensate simultaneously dis- 
appear. At small hyperradii the three-body recombina- 
tion rate is clearly much larger than the macroscopic tun- 
neling rate, whereas the opposite holds for large hyper- 
radii. For the parameters in fig. 11 we find that the two 
time-scales are roughly equal around the second mini- 
mum where the condensate is located. 

However, the tunneling rate depends strongly on the 
barrier through the combination N\a s \/bt- Varying either 
of the three quantities would then move the tunneling 
rate up or down in fig. 11. For a larger barrier the con- 
densate would only decay by direct recombination. For 
a smaller barrier macroscopic tunneling would dominate 
and the condensate would decay by "collective" recom- 
bination of all particles in a very short time interval. 

When a few particles recombine into dimcrs and leave 
the condensate, the remaining system is no longer in an 
eigenstate of the corresponding new Hamiltonian. An 
adiabatic adjustment of Hamiltonian and wave function 
could then take place. Since fewer particles and un- 
changed a s and b t means a larger barrier, the stability 
against macroscopic tunneling of the new system is there- 
fore increased. 

This stabilization by particle "emission" could also be 
the result of the recombination in the macroscopic tun- 
neling process if the time-scale for recombination at the 
relevant small distances is longer than the adiabatic ad- 
justment time. In a possible development first a number 
of particles are emitted, the adjustments follow, and a 
larger barrier appears which traps and stabilizes the part 
of the initial wave function in the second minimum. How- 
ever, now the condensate contains fewer particles. 



C. Macroscopic collapse 

These decay scenarios are open for direct experimental 
investigations since the interaction can be changed in an 
experiment by using the Zeeman splitting to tune to a 
Feshbach resonance [12, 13, 39]. An initial value of the 
scattering length (corresponding to a stable condensate 
in the second minimum) can almost instantaneously be 
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FIG. 11: Three-body recombination rate eq. (35) in units of 
the oscillator frequency v = u)/(2n), typically of the order 
of 10-100 Hz [13], as a function of hyperradius for N = 100, 
a s /b — —50, bt/b = 10 4 . Shown as the horizontal, dashed 
line is the macroscopic tunneling rate eq. (38). Shown as 
the horizontal, dotted line is the macroscopic collapse rate 
eq. (44) when the scattering length is much larger than the 
trap length. 

changed to a value where the barrier is removed. The 
initial wave function for the condensate is now no longer 
a stationary state in the new potential. 

If we assume that the only excitations are the degrees 
of freedom contained in the lowest new hyperspherical 
potential with s-waves, we can use the sudden approxi- 
mation and expand on the corresponding cigenfunctions. 
The most important levels are then the lowest-lying posi- 
tive energy states with energies comparable to the initial 
condensate. The time-scale for the time evolution of the 
initial state in the new potential is then determined by 
the energy differences between such levels. These states 
of positive energy and large spatial extension confined 
by the trap are roughly separated by the oscillator quan- 
tum of energy two. The corresponding rate for populating 
smaller distances with the consequence of immediate re- 
combination is then crudely estimated to be 



Experimentally [13] this macroscopic collapse time is ver- 
ified to be of the order ~ l/u>, typically a few millisec- 
onds, as given by the external trapping field. 

This macroscopic collapse time is shorter than the 
macroscopic tunneling time for the parameters of the sys- 
tem in fig. 1 1 . The motion in the potential is fast or slow 
compared to the recombination time for distances in the 
first or second minimum, respectively. The time evolu- 
tion after the sudden removal of the barrier could then be 
a macroscopic collapse towards smaller hyperradii where 
dimers and trimers are "emitted" and the barrier begins 
to appear. The part of the wave function trapped at 
large distances in the second minimum can then stabilize 
into a condensate with fewer particles. The time-scale 



for these processes should then be between the macro- 
scopic collapse time and the recombination time at the 
second minimum. Possibly other time scales due to the 
neglected degrees of freedom (angular momentum, clus- 
terization, etc.) could be present in the foil study of the 
dynamics of a many-boson system. 



VI. SUMMARY AND CONCLUSION 

Correlations in a system of N identical bosons of low 
density are described by use of a hyperspherical adia- 
batic expansion. The wave function is decomposed in 
additive Faddcev-Yakubovskii components, where each 
term is related to one pair of particles and only s-waves 
are included. The adiabatic potentials are only weakly 
coupled and we investigate structures where only the low- 
est contributes. We use a finite-range purely attractive or 
purely repulsive Gaussian interaction and extract general 
properties of the lowest angular eigenvalue. 

We establish universal scaling relations for the ra- 
dial potential for arbitrary scattering length and particle 
number. These scaling rules are valid for large and inter- 
mediate distances where the particles on average are out- 
side the range of the interaction. Only the short-distance 
behaviour is influenced by the choice of interaction po- 
tential. 

We parametrize the model-independent part of the ef- 
fective radial potential in a simple form with an inter- 
action part, a centrifugal barrier term and a contribu- 
tion from the external field. This potential diverges at 
small distances due to the centrifugal barrier and at large 
distances due to the confining external field. The two 
minima are generally separated by a barrier. The deep- 
est minimum at small to intermediate distances supports 
self-bound iV-body systems where the density is much 
larger than for a Bose-Einstein condensate. The second 
minimum at a much larger distance allows solutions with 
properties characteristic of a condensate. We distinguish 
by formulating a definition of a condensate in this con- 
text. 

We compare properties of the correlated structures 
with those of the zero-range mean- field solutions. The 
large-distance asymptotic behaviour is found numerically 
to reproduce the mean-field result for a zero-range inter- 
action renormalizcd to give the correct scattering length 
in the Born approximation. This is remarkable since the 
correct scattering length for the Gaussian potential is far 
from the Born approximation. Thus the different terms 
in the second-order integro-differential equation conspire 
to produce this large-distance result, which is rigorously 
established for three particles and on general grounds also 
expected for many particles. The choice of wave function 
is then a posteriori shown to be sufficient. 

The stability of the condensate is limited by decay 
into lower-lying many-body cluster states reached by 
processes where three-body recombination resulting in 
bound dimers is very prominent. We compute various 
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rates of decay and discuss the time-scales involved. The 
bare three-body recombination process is strongly scatte- 
ring length and density dependent and therefore increases 
dramatically when the wave packet moves from the se- 
cond minimum to smaller distances. An intermediate 
barrier would only allow quantum tunneling followed by 
a macroscopic collapse. When this barrier is very small 
by choice of parameters the macroscopic tunneling rate 
would dominate. When the interaction is changed dur- 
ing an experiment and the barrier is totally removed the 
already created condensate would collapse and a number 
of cluster configurations would appear. Stability may 
subsequently be automatically restored and a new con- 
densate created with fewer particles. 

In conclusion, we have discussed properties of conden- 
sates and extracted universal scaling relations. We have 
focused on the effects of correlations for large scatte- 
ring lengths where the mean-field approximation breaks 
down. Finally we investigated time-scales for various de- 
cay mechanisms limiting the stability of the condensate. 
The parametrized potentials allow independent investiga- 
tions without the full numerical machinery. More general 
TV-body structures are studied than the simple conden- 
sates. 



APPENDIX A: NUMERICAL DETAILS 

The angular equation can be scaled by using the poten- 
tial range b as the unit length [27]. The only interaction 



parameter is then the Born approximation to the scatte- 
ring length in this unit a^/b. The only length coordinate 
is then p/b. All physical quantities are functions of such 
dimcnsionlcss ratios. 



The s-wave two-body scattering length is the node of 
the zero-energy solution to the two-body Schrodingcr- 
equation, i.e., u(r) cx (r — a s ). Table II shows the scatte- 
ring length a s for different potential strengths <zb, see 
cq. (16). The Born-approximation equals the correct 
scattering length only in the limit of weak attraction, 
where the scattering length a s is much smaller than the 
range of the interaction b. 



To exemplify, in experimental work Rb atoms with 
a scattering length of a s = 100 a.u. are trapped in an 
external trap of frequency v = 100 Hz [22]. Assuming an 
interaction range around b = 1 nm we obtain a s /b = 5.29, 
b t /b = 1442. This can be modelled by a Gaussian two- 
body interaction with a-^/b = —1.5, where the lowest 
solution corresponds to a two-body bound state and the 
next accounts for the properties of the condensate. 
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TABLE II: The scattering length a s in units of b for various 
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